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This paper describes a statistical method for short-term forecasting of surface layer 
wind velocity amplitude relying on the notion of continuous cascades. Inspired by 
recent empirical findings that suggest the existence of some cascading process in the 
mesoscale range, we consider that wind speed can be described by a seasonal compo- 
nent and a fluctuating part represented by a "multifractal noise" associated with a 
random cascade. Performances of our model are tested on hourly wind speed series 
gathered at various locations in Corsica (France) and Netherlands. The obtained 
results show a systematic improvement of the prediction as compared to reference 
models like persistence or Artificial Neural Networks. 
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I. INTRODUCTION 

The fast growth of wind energy technology shows that more and more countries attach im- 
portance to this renewable resource. However, the energy production is strongly dependent 
on the wind volatility and is consequently characterized by a large amount of uncertainty. 
Reliable wind speed predictions are therefore necessary to optimize plants scheduling or to 
evaluate systems production. For that purpose, many efforts have been spent for several 
years by the scientific community in order to design faithful models that allow one to perform 
good forecasts. As reviewed e.g. in there are mainly two families of approaches. The 
"physical" models rely upon physical considerations leading to some atmospheric models 
that provide a "numerical weather prediction" system. For very short prediction horizons, 
one often prefers "statistical" approaches that mainly consist in designing stochastic models 
or using methods of time series analysis, calibrated on historical data or other explanatory 
variables (like the output of a physical model). Within this framework, one can cite standard 
ARIMA modeling j2n4|, models relying on Markov chains jsf, wavelet based methods j^, 
"black boxes" methods like advanced Recursive Least Squares or Artificial Neural Networks 
(ANN) 0. 

The method we propose in this paper is based on recent empirical results according to 
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which short time wind variations possess intermittent statistical properties similar to those 
actually observed in fully developed isotropic turbulence (sf : they are strongly non Gaussian 
and characterized by long range correlated (log-) amplitudes |9|. These features have been 
shown to be the hallmark of random cascade processes. We therefore propose to build a time 
series wind speed model involving a multifractal noise. We aim at performing predictions of 
the wind intensity over horizons extending from 1 hour to 48 hours, using various data series 
recorded at different sites located in Corsica (France) and Netherlands. We then compare 
our results to those obtained with other common forecasting methods such as persistence, 
often considered as a reference, or an Artificial Neural Network. 

The paper is structured as follows : in section |TT] are described the various time series 
used in this study. After a brief review of their main linear properties (power spectrum, 
seasonality and correlations), we recall the observations of Muzy et al. [9| concerning the 
statistics of wind variations amplitude. Section IIIII is devoted to the definition of a simple 
stochastic multifractal model for the wind velocity components relying on former observed 
features. In section IIVI we present results of the application of this model to short term 
predictions and comparison to the aforementioned reference predictors. Conclusion and 
prospects are provided in section |Vl 



II. DESCRIPTION OF THE WIND SPEED TIME SERIES 



A. Presentation of the data and basic statistical properties 



Location 


Latitude 


Longitude 


Dates 


Sampling freq. 


Site 


Vignola (Ajaccio) 


41°56'N 


8°54'E 


1998-2003 


1 min 


70m, coastal, high hills 


Ajaccio 


41°55'N 


8°47'E 






5m, coastal, plain, airport 


Bastia 


42°33'N 


9°29'E 






10m, coastal, plain, airport 


Calvi 


42°31'N 


8°47'E 


2002-2006 




57m, coastal, hills 


Conca 


41°44'N 


9°20'E 




1 hour 


225m, high hills 


Renno 


42°11'N 


8°48'E 




755m, mountains 


Sampolo 


41°56'N 


9°07'E 






850m, mountains 


Eindhoven 


51°44'N 


5°41'E 


1960-1999 




20m, plain 


Ijmuiden 


52°46'N 


4°55'E 


1956-2001 




4m, coastal, plain 


Schipol 


52°33'N 


4°74'E 


1951-2001 




-4m, plain, airport 



TABLE I: Main features of the time series 



The time series used in this paper are amplitude and direction of horizontal wind speeds 
recorded in Corsica (France) and Netherlands. The first series called "Vignola", has been 
recorded (at 10 meters height) by the means of a cup anemometer, every minutes during 
6 years (1998-2003) at our laboratory (Vignola) in Ajaccio, Corsica Island, France. Other 
data sets consist in five years horizontal wind speed, determined every hour (10 minutes 
averages) for 6 sites in Corsica. These data have been measured and collected by the french 
Meteorological Service of Climatology (Meteo-France) using a cup anemometer and wind 
vane at 10 meters above ground level. For comparison purpose, we have also studied the 



wind data freely available from KNMI HYDRA PROJECT [10| . These data represent series 



of hourly mean potential wind speeds recorded in 3 different sites in Netherlands. Table |T] 
summarizes the main features of the sites. 
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FIG. 1: Power spectrum density of wind intensity x-component series : Vignola at the top and 
Eindhoven below, in log-log representation. 

In the sequel, V{t) will denote the modulus of the velocity horizontal vector while Vx{t) 
and Vy{t) will stand for its two components along arbitrary orthogonal axes x and y. We 
have by definition: 



vit) = ^v.{ty + Vy{ty 

Since our goal is to construct a parsimonious stochastic model of wind variations, we have 
chosen to study and Vy separately, a modeling of the wind dynamics in polar coordinates 
(modulus and direction) would be cumbersome and more difficult to handle using Gaus- 
sian processes (see next section). Moreover, since there is no well defined wind direction 
with a small "turbulent rate", components along longitudinal and transverse directions are 
meaningless and we have preferred to focus on components along arbitrary fixed directions. 

The power spectrum is one of the most common tools for analyzing random processes. 
Since the pioneering work of Ven Der Hoven [ill, HI, the shape of a typical atmospheric 



wind speed spectrum in the atmospheric boundary layer is still matter of debate. It is rela- 
tively well admitted that it possesses two regimes separated by low energy valley called the 
"spectral gap" located at frequencies around few minutes. This gap separates the microscale 
regime, where turbulent motions take place, from the mesoscale range. In Fig. [T]are plotted, 
in log-log representation, the power spectrum of Vx wind component series corresponding to 
Vignola and Eindhoven sites (for the sake of clarity, the graphs have been shifted by an arbi- 
trary constant). One sees that the Vignola spectrum (top curve) allows one to resolve higher 
frequencies than the Eindhoven spectrum. In the former one, the beginning of the spectral 
gap "plateau" can be clearly observed while the Eindhoven series goes down to smaller fre- 
quencies since it covers a wider time period. The first striking feature of both spectra are 
the main peaks associated with diurnal oscillations (see below). Up to the presence of these 
peaks, both spectra can be represented by a decreasing function that connects the fiat low 
frequency behavior to the high frequency spectral gap. The exact shape of the spectrum 
in the (intermediate) mesoscale range is unknown but it can be modeled by a power-law 
P{f) ~ with an exponent /3 between 1.5 and 2. One does not expect the same level 
of universality of the speed statistics as for turbulence at mesoscale range and notably the 
value of the exponent can depend on local orographic, atmospheric conditions.... (l3|. Let 
us note that the spectrum associated with the Vy velocity component behaves in a similar 
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way. The main power spectrum features can be alternatively observed through the behavior 
of the correlation function of the velocity components. In Fig. [2] is plotted the estimated 
covariance of the de-seasonalized component of Schipol wind data as a function of the 
lag r (see section [ITIl for the details about the way we process seasonal effects). One sees 
that the correlation decreases quite slowly and the velocity remains correlated up to lags of 
few days. Wind components auto-correlations and cross-correlations can be easily described 
within the framework of linear time-series models. ARMA like modeling fl3] have been 
widely used to model many meteorological time series [l5| like monthly precipitation 16|. 
annual streamflow 17 1 or monthly drought index 18[. Many authors have also considered 
such time series models in order to account for the fluctuations of the wind velocity ampli- 
tude or its components (see e.g., [l, 0, 



19|-|22[). However, since most of these approaches 



are faced to the non Gaussian nature of the wind fluctuations (see next section) and the 
presence of seasonal effects, many of these models involve some non- linear "normalization" 
transformation and/or a separate parametrization of each season |i,[23|-[25|. It results that, 
despite the simplicity of ARMA processes, the final models remain relatively hard to esti- 
mate and far from being parsimonious. In this paper we choose to use (seasonal) ARMA 
processes and account for the non-gaussian observed statistics through the nature of the 
noise term that will be given, as explained in the next section, by a multifractal process. 
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FIG. 2: Estimated covariance of deseasonalized Vx component for Schipol (Netherlands). 



B. Non-linear statistical properties : non gaussian fluctuations and magnitude 

long-range correlations 

When referring to the non-Gaussian nature of wind speed statistics, one has to be precise 
since velocity amplitudes are obviously not normally distributed. Indeed, even in the case 
when Vx and Vy are Gaussian, the velocity modulus pdf is a Rayleigh distribution (or a Rice 
distribution if the components have a non zero mean) a particular case of the Weibull family. 
This is probably the main reason why Weibull is the most commonly considered distribution 
in order to reproduce the pdf of wind amplitudes [26|. All the approaches that consist 
in directly trying to describe the stochastic dynamics of the wind amphtude are faced to 
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FIG. 3: Probability density function (pdf) of the error px observed for an AR(1) model of the 
component of Schipol wind series. The dashed line corresponds to a standardized normal law. 



problems related to the non Gaussian nature of its statistics. In particular, when one wants to 
account for the observed linear correlations, since the involved distributions are not stable by 
aggregation (unlike Gaussian random variables), a control of both persistence (correlations) 



and the nature of the statistics is very difficult [21|, |27H29[ . As mentioned previously, some 



authors have tried to reproduce the wind correlations, within AR Gaussian models, by using 
a non-linear transformation of wind amplitudes in order to handle normal random variables 
[21, |3| . However, beyond the fact that these methods make strong assumptions on the nature 
of empirical laws, they only account for the mono-variate distribution ^ and they are not 
stable as respect to aggregation, in the sense that a change of the sampling period or the 
size of time averages would drastically modify the parameters involved in the model. 

As discussed in the previous section, a modeling of both components Vx and Vy allows one 
to reproduce the observed (partial) correlation functions within the framework of ARMA 
models. However, even in the context, non Gaussian statistics are observed: if one studies 
the distribution of the prediction errors of these models (or simply the distribution of velocity 
components variations), it appears that their pdf are characterized by stretched exponential 
tails very similar to the distribution of velocity increments at small scales in fully developed 
turbulent flows jsf. In Fig. [3] is plotted the logarithm of the standardized pdf of the noise 
obtained using an AR(1) process in order to model the variations of the Schipol series Vx 
component (the additive seasonal part has been removed). For comparison purpose we have 
also plotted the parabola associated with a standardized Gaussian law. It appears clearly 
that the pdf of noise fluctuations has a kurtosis very large as compared to the normal law. 

Another striking feature of wind series is that the amplitude of the error noise is long-range 



correlated. This is another property that has been observed in turbulence |30l . |31| . More 



t Indeed, even if marginal probability densities are Gaussian, nothing guarantees that it is the case for the 
n-variate distributions 
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FIG. 4: Square root of wind magnitude covariance as a function of the log of the lag (time units 
are hours). The solid line corresponds to data from Vignola (20 minutes average). The symbol (□) 
curve represents the average over the 7 Corsica sites, and the symbol (■) curve, the mean of the 
3 sites from Netherlands. 

precisely, if one defines the local 'magnitude' as z/(t) = |ln(p^.(t)^) or z/(t) = |ln(pj^(t)^) or 

uit) = hn{p,ity + py{ty), (1) 

where px and py are the noise terms associated with a linear prediction of Vx and Vy ^, then 
the empirical covariance of u can be fitted as: 

CovHt),u{t + r)]^^Hn'(^^) . (2) 

By representing the square root covariance as a function of the logarithm of the lag r one 
should obtain a straight line. This is illustrated in Fig. HI where the magnitude covariances 
estimated for the sites in Corsica and Netherlands have been plotted (see caption). It can 
be observed that the parameters (slope of the curves) and T (time lag where correlations 
vanish) are very close for all site data. As explained in Ref. ^ and briefly reviewed in the 
next section, these properties are intimately related to random cascade processes. 



C. Random cascade model for vi^ind speeds 

Discrete random multiplicative cascades were originally introduced as models for the 
energy cascade in fully developed turbulence. In the simplest case, these objects are positive 



^ We would obtain the same results with i/s(i) = \n{SsV^ + SsVy)/2 where s is a small scale and SsF{t) = 
F{t + s) - F{t) 
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fields (measures) whose construction involves a recursive procedure along a dyadic tree: the 
cascading process starts at a large "integral" time scale T where the measure is uniformly 
spread (meaning that the density is constant). One then splits this interval in two equal parts 
over which the densities are obtained by multiplying the 'father' density by two (positive) 
i.i.d. random factors Wi = e'^^ and W2 = e"^. Each of these two sub-intervals is again cut 
in two equal parts and the process is repeated infinitely. At construction step n, the dyadic 
intervals have a size T2~"' and their measure denoted cr^ is simply: o"^ = o"o 11 ^« ~ 2~"e^''% 
where all the Wk = e'^* are i.i.d such that E{W) = 1. If the random variables k are 
Gaussian, then the corresponding model is log-normal and its scaling properties are easy 
to control (see e.g., [s^ and references therein for more details). Let us notice that non 
positive fields like, for example, the velocity field in developed turbulence, can be simply 
derived from the construction of the measure by considering that o"^(t) is the (stochastic) 
variance of a Brownian motion (or another Gaussian process), i.e., 5T-X{t) = ^/a^{T)e , 
where e is a Gaussian random noise. Such "grid bounded" cascades, though simple, do not 
however provide a satisfying model for a stationary physical process such as wind temporal 
variations. Indeed, they are built on a fixed time interval [0,T], are not causal and not 
stationary. Moreover, they involve an arbitrary fixed scale ratio (2 in the dyadic case). 
Very recently, several constructions have been proposed to generalize discrete cascades to 



stationary, causal and continuous processes [32|, |33|. We will not enter into details but if 
one calls the magnitude process uit) = ^ ■ Kj, then in the log-normal case, w is a gaussian 
process characterized by its covariance function. If one notices that the tree-like structure 
underlying the discrete construction implies a logarithmic correlation function, then one can 



naturally define the log-normal continuous cascade as follows |33L l34j| : 

al{t) = e^'^^W (3) 
where Ws(t) is a stationary gaussian process of covariance defined by : 

CoY[us{t),Us{t + T)] = XHn{^). (4) 

S + T 

Here T and are two parameters that correspond respectively to the integral scale (corre- 
lation length analog to the time scale where cascading process starts) and the intermittency 
coefficient (which quantifies the degree of burst occurrences in the process). The parameter 
s is a time sampling parameter that can be chosen arbitrary small (since the weak limit 
s — )■ of the process exists j32|, [1^1). It can be proven that such a process is the continuous 
equivalent of discrete random cascades. Therefore, according to this picture, a continuous 
cascade is nothing but a stochastic process which magnitude, as defined by the logarithm of 
its variations, has a covariance correlated as a logarithmic function. 

In order to link these considerations with previous observed features for wind data, let us 
remark that wind fiuctuations at a fixed spatial location result from two types of stochastic 
variations: first, the spatial fiuctuations at a fixed time (Eulerian) and then the temporal 
fiuctuations for a fixed fiuid element (Lagrangian). Since there is no strong mean velocity 
and Taylor frozen hypothesis cannot be invoked, both Lagrangian and Eulerian variations 



have to be taken into account. In ref. [36j, B. Castaing shows that if one supposes a 



continuous cascade paradigm (Eq. (j3])) for both Eulerian and Lagrangian fields, then the 
magnitude correlation function at a fixed location should behave like a squared logarithmic 
function: 

Coy [uj,{t),u,{t + t)] = pHn'i-^) (5) 

S + T 
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where the coefficient depends on both Lagrangian and Eulerian intermittency coefficients. 
This is precisely the behavior that we observed in real data as reported in Fig. HJof previous 
section (see for more details). Therefore, the residual variance of errors Px{t) and Py{t) 
associated with linear models of Vx(t) and Vy(t) can be both defined as in Eq. {Px{t) = 
e'^=W£^(t) and p^(t) = e'^'^^^hy{t)) and: 

2z/(t) = ln(p,(t)=' + py{t)^) = 2uj{t) + In Z{t) (6) 
where Z{t) = e^itf + Syitf ■ 



III. BUILDING THE MODEL 

Let us now sum up all the reported empirical observations in order to build a time series 
model of wind speed components Vxit) and Vy{t). According to previous considerations, the 
model will be formulated as a seasonal auto-regressive process where errors are given by a 
(seasonal) continuous cascade. 

A. Construction of the seasonal autoregressive part 

It has been shown in section fll Al that Vx{t) and Vy{t) both contain additive seasonal 
components, i.e., can be written as: 

Vx,y{t) = Sx,y{t) + Vly{t) (7) 

where Sx,y{t) represent the deterministic diurnal oscillations and V^y{t) the "de- 
seasonalized" velocity components. Since the seasonality is caused by the variation of the 
sun position during the day, Sx^yif) are almost daily periodic functions, with a period shape 
that changes according to the considered season in the year. In order to determine this 
shape, we therefore have to perform a "local" estimation. For that purpose, we use a stan- 
dard methodology described in [37]: each seasonal component Sxif) and Sy{t) (denoted as 
S{t)) is described by m Fourier modes of period 1 day {D = 24 samples for hourly data): 

m 

S{t) = ao + Y, 

k=l 

Because of the yearly variation of the seasonality, the coefficients {?7i,fc}j=i,2;fc=i...m depend a 
priori on the day d and the local estimation simply consists in using least squared method 
associated with a local exponential moving average: 

{Y D-l 
^^^M-.l^[V;,,(yy,j,t)-5(t)f 
yy=l j t=0 

where Y is the number of available years in the data series, Vx^y{yy,j, t) represent the velocity 
component at year yy, day j and 'hour' t. is an exponential discount factor chosen so 
that j^j^ ~ 10 days {ip = 0.9). We have used D = 24 for hourly data and m = 3. We have 
checked that our results remain almost unchanged if one increases the number of harmonics. 



.2kTTt. .2k7it, 
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FIG. 5: PACF versus time lag : (a) corresponds to Schipol data and (b) to Ajaccio data. 

Empirically, we have found that seasonal components represent 20 to 35 % of the wind 
amplitude energy, except for 2 sites, Ajaccio and Renno (Corsica) where they represent 
around 50 % of the total energy. 

In order to account for the linear correlations and cross-correlations of the stationary 
parts V^(t) and Vy(t), we have considered the class of bi-variate ARMA processes. The 
study of partial autocorrelation (PACF) and cross-correlation functions suggests that an AR 
of order 2 or 3 is appropriate to fit the observations. This is illustrated in Fig. 13 where plots 
of PACF versus the lag are reported for wind speed component of Schipol and Ajaccio 
series. We have consistently observed that for all series, the PACFs are close to zero value 
after lag 2. An AR(3) model should be more appropriate for some sites, but accounting to 
higher order auto-regressive processes does not lead to any significant improvement of the 
results reported below. 

The results of correlograms study can also be confirmed by other model selection proce- 
dures like the Akaike Information Criterion (AIC). The best choice of the AR order p is the 
value that minimizes the following quantity : 

AICip) = N\n{al{p)) + 2P, (9) 

where N is the length of each data series, P is the number of estimated parameters and 
ap is the variance of the residuals p. We have studied this criterion for different sites and 
observed that AIC{p) decreases fast from p = 1 to p = 2 and slower at lags beyond 2; that 
confirms our previous results on the PACF and the choice of an AR(2) model. Considering 
orders greater than 2 does not improve the model forecasting performances. 

Finally, we are lead to the following simple model for deseasonalized wind components: 

v,'it + 1) = eLo (7.x(fc)V;^(t -k) + 7.,(fc)V;^(t - k)) + p.(t + 1) 

(10) 

Vy'it + 1) = ELo (7,,(fc)Kf ly.ik)V,Ht - k)) + Py{t + 1) 

where Px,y{t) represent the noise terms which will be modeled as a log-normal continuous 
cascade, (see next section), ■jxxik), lyyik), Ixyik) and 'jyxik) {k = 1, 2) are the AR coefficients. 

Let us notice that the values of these coefficients strongly depend on the (arbitrary) choice 
of the reference direction defining ^ and one cannot expect any universality or physical 
meaning in the precise value of each coefficient. For instance, the coefficients estimated for 
the Schipol series are 7xx(o) = 0.87, 'jxxii) = 0.09, 'jxyio) = —0.05 and 7xy(i) = 0.04 while the 
values we found for Ajaccio are 'jxx(o) = 0.56, ■jxxii) = 0.11, 7xy(o) = 0.06 and 7x1/(1) = —0.04. 
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From a methodological point of view, we split the data in two parts : the first part 
of each database (4 years for all Corsica sites, 20 for Ijmuiden, 30 for Eindhoven, 40 for 
Schipol) was used as the "training period" or the "learning part". All the parameters of our 
model, are determined using the learning part. The remaining data allow us to evaluate the 
performance of each model as it will be seen later. 

B. Accounting for the cascade 

As explained in section III CI within the random cascade paradigm, the noise p{t) can be 
written as : 

where ex,y(t) are independent white Gaussian noises, Ms(t) is a deterministic function that 
represents a multiplicative seasonality of the noise amplitude and u{t) is a zero mean sta- 
tionary gaussian sequence independent of e{t), which covariance is a squared log as de- 
scribed in section III CI (Eq. ([5])). In practice, one computes 2z/()f:) = ln(p^(t) + Py(t)) = 
2u}{t)+2Ms{t) + ln{el{t)+ey{t)) and since the mean and variance of InZ(t) = ln{el{t)+ey{t)) 
are known, one can obtain Ms{t) along the same line as we have estimated Sx,y{t) (Eq. (jHD). 
A generalized method of moments [3^ applied to the sample covariance of v{t) allows us to 
evaluate the parameters and T of Eq. ([S]), defining the Gaussian process uj{t). 

IV. APPLICATION TO SHORT TERM PREDICTION 

A. H-step forward prediction 

Our goal is to predict wind speed intensity V{t) = \JVx{tY + Vy{tY at different horizons 
of time (from 1 hour to 48 hours). Since the (conditional) law of the velocity modulus is 
not Gaussian, the "best" prediction depends, in general, on the type of error one wants 
to minimize. In theory, since the multifractal AR model we have introduced provides the 
full conditional law of each velocity component, one should be able to optimally solve any 
forecasting problem. For the sake of simplicity, we will only estimate the conditional mean 
of V , denoted as E{y\t) in the sequel, that is the predictor which minimizes the mean square 
error. 

Let V^^yit, h) (resp. V^^yit, h)) be the best linear predictors of V^y{t + h) (resp. V^^yit, h)), 
at time t and horizon h, i.e., from the definition of the model: 

Vlit,h) = E[Vlit + h)\t] 
%,y{t,h) = Vly{t,h) + Sx,y{t + h) 

These predictors are easy to compute: since the linear part of our model reduces to a vector 
AR(2) model (Eq. f lTO|) ). h iterations of the model provide the linear coefficients. Indeed, 
Eq. f llOp can be rewritten in a vector form: 



V^(t + 1) =AY^{t) + e{t+l) 



(12) 
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where the vectors V"^(t) and e(t) are defined by: 



V^(t) 



via - 1) 



, e(t) 







and the matrix A reads: 



A 



'~1xx (0) 7xy(0) 7xx(l) 7xj/(l) \ 

7yx(0) 7ra(0) 7?/x(l) 7ra(l) 
10 
10 



When one considers an horizon the iteration of Eq. f|T2|) gives: 



/i-i 



V^(t + /i) = A^\^{t) + ^A^e{t + h-k) = A^'Y^it) + e'^^\t + h). 



k=0 



(13) 



(14) 



(15) 



According to this representation, V^y{t,h) correspond to the first two components of 
A'^\^{t). From Eqs. (fTTl) and (fT3|) . the components of the noise vector, in the r.h.s. 
of previous equation, can be written as: 



+ ^) = E «^e^^*+^"^'^e.,,(t + h-k), 



(16) 



where the constants can be deduced from the A coefficients. Moreover, by considering, 
as shown in ref. [38|, e{t)e^^^^ quasi-stable as respect to hnear combinations, we have: 



(17) 



where e^'*-' is a standardized Gaussian noise and fi*^'^-' is also Gaussian, at fixed h, with the 
same covariance as fl{t) for lags greater than h (Eq. (|5])). Eqs. ( IT5|) and (fT7|) show that the 
model conserves the same shape for all prediction horizons: 



Vl{t + h) = Vl{t,h) + 



,nW{t+h)(h) 



iit + h). 



(18) 



This property is of great practical interest because, whatever the horizon h, at fixed value of 
^l^'^\t + h), the law of the velocity modulus V{t + h) is a Rice distribution 39| of parameters 

r = ^Jv^{t + h) + Vy^{t + h) and = e^^'^'^'' More specifically, let M^(r,a2) be the 
mean value of a Rice distribution, i.e., 

.2 



r 



2a2 



(19) 



(where Li/2{x) is the order 1/2 Laguerre polynomial), and Ph{^\t) the conditional Gaussian 
law of il^^\t + h). The conditional velocity value at horizon h is then: 



E{V{t + h)\t) = I Ph{VL\t)MR{r,e^^\t)dVL 



(20) 
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This quantity can be evaluated numerically by a Gaussian quadrature approximation of the 
Gaussian integral 40|. The conditional law of ^l^^\t + h) is a normal law which mean, 
Q^''\t + h), and variance, SQ\t-\-h), can be computed using the known mean and covariance 
of Q^^\ Cl^^\t + h) is nothing but the best linear predictor of + h) at time t and horizon 
/i, i.e.: 

T-l 

QS''\t + h) = MP{t + h) + ^akJ''\t-k), (21) 

fc=0 

where the filter size T and the coefficients are obtained from the shape of the covariance 
function of uo^^'^ (Eq. (ED). If one denotes C\f = Gov [u^''\t) , u'^^\t + \j - and Cf ^ = 
Gov [u'^^\t),u^''\t + k + h)], then 



EK'fcf- (22) 

j 

Let us end this section by noticing that the alternative predictor 



Vit + h) = ^/E{V^{t + h)\t) , (23) 
which, after a little algebra, reduces to 



V{t + h) = \lv^{t + hf + Vy{t + hy + 2e2f^^''^(*+'^)+24"(*+'^) (24) 
provides performances relatively close to the former "Rice" predictor. 



B. Forecasting performances of our wind model 

We present in this section the forecasting performances of the previously defined model 
as compared to standard models like persistence, a reference model introduced by Nielsen et 
al. [llj and a simple Artificial Neural Network (ANN). The parameters of these two latter 
models are estimated over the previously defined "learning part" of each data series (see 
section iniAl). Models comparison are made using two different mean error measurements. 



1. Reference models 

As explained in [l| , simple techniques are often used as references within the wind power 
forecasting community. Let us briefly describe the 3 main models we considered for perfor- 
mance comparison purpose. 

• Persistence 

This model is the most commonly used reference predictor. According to Giebel [42|, for 
short prediction horizons (from few minutes to hours), this model is the benchmark all other 
prediction models have to beat. It consists in a simple martingale hypothesis according to 
which future wind speed at horizon h will be the same as the present observed value : 



V{t + h\t) = V{t). 
Merge of persistence and global average 



(25) 
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Nielsen et al. 4l|] propose to use a linear combination of the persistence predictor and 
the global average to improve the previous persistence prediction: 



V{t + h\t) = aV{t) + il-a)V, 



(26) 



where a is the correlation coefficient between V{t) and V{t + h) and V is the mean velocity. 
V and a can be determined using data up to time t or using the chosen training period of 
each database. Let us note that V — V can be identified as an AR(1) predictor. 
• Artificial Neural Network (ANN) model 

Artifical neural networks are commonly used as "black boxes" prediction tools in many 
areas. Notably, there is a wide literature on their interest in wind speed forecasting (see 
e.g. [l| and references therein). We have designed this method using the ANN toolbox of 
MATLAB, with the collaboration of Philippe Lauret, as introduced in [43[. We have chosen 
the most popular form of NN called multilayer perceptron (MLP) structure. The MLP 
structure consists of an input layer, one or several hidden layers and an output layer. In our 
case, the input vector is given by the previous observed values of the wind speed and the 
output vector consists of only one output, which is the corresponding forecast at horizon h. 
Best results are here obtained with 30 input neurons and one hidden layer, characterized by 
5 non-linear units (or neurons). The non-linear function associated with each unit is usually 
a tangent hyperbolic function f{x) = tanh(x). Therefore, a NN with Ni = 30 inputs, Nh = 5 
hidden neurons and a single linear output unit defines a non-linear parameterized mapping 
from an input x to an output y, given by : 



y = y{x;w) = 



j=0 



Wji.Xi 



(27) 



where wj are the weight applied on each hidden neuron and Wji ones applied on each input 
data. These NN parameters w are estimated during a learning phase. It consists in adjusting 
w so as to minimize an error function which is usually the sum of squares error between 
measured data and network output (see next section). For that purpose, several iterations 
are necessary (we have observed that 30 are sufficient). 



2. Estimation of forecasting accuracy 

Errors frequently used to compare various prediction methods are the Root Mean Square 
Error (RMSE), the Mean Absolute Error (MAE), the Mean Error (ME), histograms of the 
frequency distribution of the error or the correlation function We have chosen to employ 
the most common of them, i.e., the RMSE and the MAE, in order to evaluate the relative 
performances of each model. These errors are given as percent of the mean of wind speed at 
each site. If V(t) is the observed wind speed at time t and V{t) the corresponding forecast, 
these errors are defined as follows: 

• The normalized mean absolute error (nMAE) is simply defined as: 

1 1 " 

nMAE = =-J2\ " I' (2^) 

^ ^ t=i 

where n is the number of periods of time and V is the mean velocity amplitude over the 
testing period. 
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• The normalized root mean square error (nRMSE), which gives more weight to largest 
errors, reads: 

nRMSE = LvMSE (29) 



V 

with 



1 

MSB = - - V{t) f (30) 



n 
t=i 



3. Results 



Location 


Pers. 


Pers+Mean 


RNA 


Mult. mod. 


Vignola 


42.7 


39.8 


38.4 


37.6 


Ajaccio 


40.4 


36.6 


34.8 


33.8 


Bastia 


44.9 


42.1 


40.4 


40.2 


Calvi 


40.2 


38.4 


35.7 


36.0 


Conca 


49.7 


47.4 


46.0 


46.0 


Renno 


44.1 


40.5 


39.1 


37.6 


Sampolo 


54.4 


51.6 


48.3 


47.9 


Ijmuidou 


13.G 


13.5 


13.5 


13.6 


Schipol 


17.5 


17.3 


17.1 


16.9 


Eindhoven 


20.3 


20.0 


19.8 


19.7 



TABLE II: nRMSE (%) of each site at one hour horizon. Best results are indicated using bold 
faces. 



Location 


Pers. 


Pers+Mean 


RNA 


Mult. mod. 


Vignola 


70.3 


56.2 


53.5 


51.2 


Ajaccio 


66.6 


48 


43.1 


41.4 


Bastia 


77.1 


61.8 


57.5 


55.3 


Calvi 


66.2 


57.7 


54.4 


52.0 


Conca 


78.7 


69.1 


66.7 


66.4 


Renno 


71.3 


54.9 


52.4 


49.6 


Sampolo 


101.9 


81.8 


69.4 


65.4 


Ijmuiden 


33.5 


31.6 


31.4 


31.3 


Schipol 


43.6 


40.2 


38.7 


36.7 


Eindhoven 


47.6 


43.5 


41.5 


39.5 



TABLE III: nRMSE (%) of each site at 6 hours horizon. Best results are indicated using bold 
faces. 
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10 20 30 40 10 20 30 40 

horizon (hours) horizon (hours) 



FIG. 6: Evolution of RMSE and MAE values for different models depending on the horizon (1 
hour to 48 hours). Figures (a) and (b) illustrate these evolutions for Ajaccio (Corsica), (c) and 
(d) correspond to Eindhoven in Netherlands. For each case, symbol (□) curve represents the 
persistence model's results, (■) curve, the merge of persistence and global average, (o) curve, the 
ANN model and (•) curve, the multifractal model. 

In Fig. |6]the nMAE and nRMSE associated with each model prediction are represented, 
at various forecasting horizons, for 2 sets of data (Ajaccio in Corsica and Eindhoven in 
Netherlands). For the 1 hour horizon, as it can also be observed for the nRMSE in table 
IIVB31 the performances obtained with the cascade model are slightly better than those 
obtained with concurrent models (average improvement of respectively 1 and 10 percent as 
compared to ANN and persistence). When the horizon increases, the performances of each 
model decrease, but the relative accuracy of our model becomes more and more significant. 
This is confirmed in table IIVB3I where are reported the nRMSE at 6 hours horizon for all 
the data series (average improvement of respectively 4 and 26 percent as compared to ANN 
and persistence). 

We have also evaluated the models performances when one increases the sampling fre- 
quency of the data used to compute the prediction, for some fixed time scale and horizon. 
The data set gathered at Vignola, sampled at 1 minute rate, allows us to compare the 
forecasts of hourly mean velocity, 1 and 6 hours ahead, by using velocity data at different 
sampling rates: 10 minutes, 20 minutes, 30 minutes and 1 hour. In Fig. [7] are reported the 
prediction errors of the cascade model as a function of the sampling rate for fixed horizon 
and averaging time scale. One clearly sees a systematic improvement of the accuracy as one 
uses better resolved input data: the finer the sampling rate, the better the forecast. Similar 
improvements can be observed with others models. This result highlights the importance of 
having high frequency data to enhance the forecast quality. 
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20 40 60 

Sampling rate (nun) 



FIG. 7: nRMSE values evolution of the hourly mean Vignola series forecast using multifractal 
model, depending on the database sampling rate (10 minutes to one hour). Figure (a) illustrates 
this evolution for one hour horizon forecast whereas (b) corresponds to 6 hours horizon forecast. 

V. CONCLUSION 

In this paper, we have addressed the problem of short term wind speed forecasting using 
a simple autoregressive seasonal model involving multifractal fluctuations. This model re- 
lies on "universal" empirical observations showing that high frequency velocity components 
variations have long range correlated amplitudes [9j]. Our model is relatively parsimonious 
and accounts for the wind properties over all time scales. It has been applied to forecast 
hourly wind speed data up to two days (48h) ahead. The obtained results show that the 
proposed method is more accurate than standard reference models. Let us notice that our 
approach can be improved by considering, for instance, its natural multivariate generaliza- 
tion. This may allow us to describe the joint wind variations at different locations. Let us 
also mention that unlike 'black boxes' approaches, our time series cascade model is able to 
provide unconditional and conditional velocity probability distributions and therefore ad- 
dress many questions related to resource assessment or risk management. This problem will 
be the scope of a further study. 
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